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ABSTRACT 

The  technique  developed  by  Arons  (1948)  for  determining  the  parameters 
that  characterise  the  adiabatic  model  of  underwater  explosion  bubble  contents 
is  revised  to  account  for  the  effect  of  non-spherical  collapse.  Additional  modi¬ 
fications  are  made  to  enhance  accuracy.  The  parameters  predicted  are  similar 
to  those  given  by  Arons  but  consideration  of  the  pressure  in  the  fluid  when 
the  bubble  rebounds  shows  the  agreement  to  be  somewhat  fortuitous.  Error 
bounds  for  the  parameters  are  given  for  the  first  time. 
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The  Effect  of  Non-spherical  Collapse  on  Determination  of 
Explosion  Bubble  Parameters 


EXECUTIVE  SUMMARY 

Over  recent  years  there  has  been  considerable  interest  in  the  behaviour  of  the  bubble 
formed  when  a  mass  of  explosive  is  detonated  underwater.  Of  specific  concern  is  the 
behaviour  of  the  bubble  upon  collapse,  when  a  high  speed  liquid  jet  is  formed,  generally 
in  a  direction  towards  nearby  structures.  Although  it  is  possible  to  predict  the  shape 
of  bubbles  as  they  collapse,  no  decisive  computational  capability  has  been  demonstrated 
for  predicting  the  loads  experienced  by  a  target  due  to  nearby  explosion  bubble  collapse. 
Attempts  are  underway  to  utilise  boundary  integral  methods  to  predict  target  loadings. 

A  key  challenge  in  application  of  the  boundary  integral  method  is  to  determine  the 
parameters  that  describe  the  gaseous  bubble  contents  in  terms  of  the  mass  and  type 
of  explosive,  and  the  depth  at  which  it  is  detonated.  It  is  typically  assumed  that  a 
quantity  of  ideal  gas  is  contained  within  the  bubble  and  undergoes  adiabatic  expansions 
and  contractions  as  the  bubble  oscillates.  Such  a  model  is  characterised  by  the  bubble 
energy  per  mass  of  explosive,  denoted  by  e,  the  ratio  of  specific  heats,  7,  and  a  constant, 
/«,  that  defines  the  adiabatic  along  which  the  gas  expands  and  contacts. 

Methodologies  for  determining  these  three  parameters  are  generally  based  upon  mea¬ 
surements  of  key  physical  quantities  associated  with  the  bubble  oscillation,  with  the  con¬ 
tribution  of  Arons1  essential  in  this  regard.  Arons  embarked  upon  a  determination  of  e, 
7  and  k  by  assuming  the  bubble  remains  spherical  throughout  the  first  oscillation  period. 
This  assumption  allows  analytical  expressions  to  be  obtained  for  the  period  and  peak  pres¬ 
sure  caused  by  the  first  bubble  rebound.  The  three  parameters  are  determined  so  that 
spherical  bubble  model  predictions  are  in  close  agreement  with  measurements  of  period 
and  peak  pressure  at  rebound. 

Investigations  of  recent  years  have  made  it  clear  that  the  assumption  the  bubble  re¬ 
mains  spherical  throughout  its  motion  is  adequate  during  the  growth  and  early  collapse 
phases,  but  fails  as  the  bubble  approaches  minimum  volume  and  rebounds.  Experimental 
and  theoretical  studies  have  shown  the  formation  of  a  high  speed  liquid  jet  upon  collapse, 
which  travels  through  the  bubble  and  causes  it  to  evolve  into  toroidal  form.  For  deep 
detonations  the  jet  is  directed  upwards  owing  to  the  action  of  buoyancy. 

The  principal  contribution  of  this  report  is  to  modify  the  technique  of  Arons  to  ac¬ 
count  for  the  effect  that  deformation  from  spherical  shape  has  upon  the  pressure  field  and 
associated  determination  of  e,  7  and  «.  Error  bounds  for  these  parameters  are  also  given 
for  the  first  time.  Having  determined  e,  7  and  «,  the  variation  of  peak  pressure  at  bubble 
minimum  with  charge  mass,  distance  from  the  charge  and  depth  is  investigated.  The  be¬ 
haviour  is  considerably  different  from  that  predicted  using  the  spherical  bubble  model  and 
is  interpreted  in  terms  of  the  non-spherical  collapse  behaviour  of  the  bubble.  In  view  of 
the  complex  behaviour  shown,  it  appears  to  be  somewhat  fortuitous  that  Arons  obtained 
good  estimates  of  the  three  parameters  on  the  basis  of  the  spherical  bubble  model. 

lJ.  Acoust.  Soc.  Am.,  20,  277-282,  1948. 
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1  Introduction 

Over  recent  years  there  has  been  considerable  interest  in  the  behaviour  of  the  bubble 
formed  when  a  mass  of  explosive  is  detonated  underwater.  Of  specific  concern  is  the 
behaviour  of  the  bubble  upon  collapse,  when  a  high  speed  liquid  jet  is  formed,  generally 
in  a  direction  towards  nearby  structures.  Although  it  is  possible  to  predict  the  shape 
of  bubbles  as  they  collapse,  no  decisive  computational  capability  has  been  demonstrated 
for  predicting  the  loads  experienced  by  a  target  due  to  nearby  explosion  bubble  collapse. 
Attempts  are  underway,  however,  to  utilise  boundary  integral  methods  to  predict  target 
loadings. 

This  report  concerns  the  various  parameters  that  characterise  the  boundary  integral 
calculation,  and  these  are  introduced  in  section  3.  The  key  challenge  is  to  determine 
these  parameters  in  terms  of  the  mass  and  type  of  explosive,  and  the  depth  at  which  it 
is  detonated.  The  problem  revolves  around  the  mathematical  description  of  the  gaseous 
bubble  contents.  It  is  typically  assumed  that  a  quantity  of  ideal  gas  is  contained  within 
the  bubble  and  undergoes  adiabatic  expansions  and  contractions  as  the  bubble  oscillates. 
Such  a  model  is  characterised  by  the  bubble  energy  per  mass  of  explosive,  denoted  by  e, 
the  ratio  of  specific  heats,  7,  and  the  constant  that  defines  the  adiabatic  along  which  the 
gas  expands  and  contacts.  This  constant  k  is  such  that 

PP~7  =  *,  (1) 

where  p  and  p  are  respectively  the  explosion  bubble  gas  pressure  and  density. 

Methodologies  for  determining  these  three  parameters  are  generally  based  upon  mea¬ 
surements  of  key  physical  quantities  associated  with  the  bubble  oscillation.  In  particular, 
experimental  data  for  the  first  bubble  period,  T,  and  maximum  bubble  radius,  Rm,  are 
commonly  available  for  the  explosive  TNT.  Data  for  the  peak  pressure  at  some  point  in 
the  flow  field  as  the  bubble  rebounds  for  the  first  time  are  also  commonly  available.  Arons, 
Slifko  &  Carter  (1948)  presented  such  data  for  TNT. 

In  the  context  of  this  report  the  contribution  of  Arons  (1948)  is  essential.  Arons 
embarked  upon  a  determination  of  e,  7  and  k  by  assuming  the  bubble  remains  spherical 
throughout  the  first  oscillation  period.  This  assumption  allows  analytical  expressions  to 
be  obtained  for  the  period  and  peak  pressure  caused  by  the  first  bubble  rebound.  The 
three  parameters  are  determined  so  that  spherical  bubble  model  predictions  are  in  close 
agreement  with  measurements  of  period  and  peak  pressure  at  rebound. 

Investigations  of  recent  years  have  made  clear  that  the  assumption  the  bubble  remains 
spherical  throughout  the  motion  is  adequate  during  the  growth  and  early  collapse  phase, 
but  fails  as  the  bubble  approaches  minimum  volume  and  rebounds.  Experimental  and 
theoretical  studies  (Benjamin  &  Ellis  1966,  Lauterbom  &  Bolle  1975,  Tomita  &  Shima 
1986,  Blake  &  Gibson  1987,  Vogel,  Lauterborn  &  Timm  1989,  Plesset  &  Chapman  1971, 
Guerri,  Lucca  &  Prosperetti  1981,  Blake,  Taib  &  Doherty  1986,  Chahine  &  Perdue  1988, 
Wilkerson  1989,  Best  &  Kucera  1992,  Best  1993,  Szymczak,  Rogers,  Solomon  &  Berger 
1993)  have  shown  the  formation  of  a  high  speed  liquid  jet  upon  collapse,  which  travels 
through  the  bubble  and  causes  it  to  evolve  into  toroidal  form.  For  deep  detonations  the  jet 
is  directed  upwards  owing  to  the  action  of  buoyancy.  It  has  also  been  shown  that  under 
some  circumstances  the  jet  does  not  completely  penetrate  the  bubble  before  minimum 
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volume  is  achieved  (Best  &  Kucera  1992).  In  these  cases  the  jet  continues  to  travel 
through  the  bubble  as  it  re-expands. 

The  principle  contribution  of  this  report  is  to  modify  the  technique  of  Arons  (1948)  to 
account  for  the  effect  that  deformation  from  spherical  shape  has  upon  the  pressure  field  and 
associated  determination  of  the  three  parameters  7,  e  and  k  that  characterise  the  bubble 
contents.  The  technique  is  further  modified  so  that  the  parameters  yield  predictions  of 
bubble  period  that  differ  from  the  experimental  data  over  the  range  of  their  validity  by 
the  minimum  amount,  in  a  least  squares  sense. 

In  section  2  the  experimental  data  essential  to  this  report  are  reviewed.  The  boundary 
integral  method  for  bubbles  is  reviewed  in  section  3  and  in  section  4  the  spherical  bubble 
model  is  presented  and  the  non-dimensional  parameters  that  characterise  the  boundary 
integral  calculation  are  derived  in  terms  of  7,  e  and  k.  In  section  5  the  methodology  for 
determining  7,  e  and  k  using  a  boundary  integral  method  is  described.  By  imposing  that 
the  theoretical  bubble  period  differs  from  experimental  data  by  a  minimal  amount,  in  a 
least  squares  sense,  the  parameters  7  and  e  are  obtained.  The  pressure  data  of  Arons 
et  al.  (1948)  are  then  used  to  fix  k.  As  k  is  varied  the  boundary  integral  method  is 
used  to  compute  the  pressure  in  the  flow  at  the  point  corresponding  to  Arons  et  ai.’s 
measurements,  k  is  varied  until  agreement  is  achieved  with  the  peak  pressure  at  bubble 
minimum.  Error  bounds  on  7,  e  and  k  are  computed. 

Having  determined  7,  e  and  «,  the  variation  of  peak  pressure  at  bubble  minimum 
with  charge  mass,  distance  from  the  charge  and  depth  is  investigated.  The  behaviour 
is  considerably  different  from  that  predicted  using  the  spherical  bubble  model  and  is 
interpreted  in  terms  of  the  non-spherical  collapse  behaviour  of  the  bubble.  In  view  of  the 
complex  behaviour  shown,  it  appears  to  be  somewhat  fortuitous  that  Arons  obtained  good 
estimates  of  7,  e  and  particularly  k  on  the  basis  of  the  spherical  bubble  model. 


2  Summary  of  essential  experimental  data 


Data  for  the  first  period  of  oscillation  and  maximum  radius  of  an  explosion  bubble 
have  been  compiled  over  a  range  of  charge  masses  and  depths  of  detonation.  These  data 
are  summarized  according  to  the  following  formulae: 


T  = 


K 


w1/ 3 

(£  +  &)5/6’ 


(2) 


(3) 


In  these  expressions  w  is  the  charge  mass  in  kilograms,  £  is  the  depth  of  the  detonation  in 
metres  and  £0  is  the  atmospheric  pressure  in  metres  of  water,  also  known  as  the  atmospheric 
head.  The  constants  K  and  J  are  given  by  Swift  &  Decius  (1947)  as 


J  =  3.36  m4/3kg-1/3,  K  =  2.11sm5/6kg~1/3. 


2 


It  is  pertinent  to  note  that  Swisdak  (1978)  quotes  an  identical  value  of  K  but  utilises  a 
value  of  J  =  3.50m4/3kg-1/3.  At  the  time  these  data  were  compiled  charge  mass  was 
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typically  measured  in  poimds,  and  distances  given  in  feet.  In  these  units  the  coefficients 
are  (Swift  &  Decius  1947) 

J  =  12.6  ffc^nr1/*  K  =  4.36  s  ft5/6lb-i/3) 

with  Swisdak  giving  J  =  13.1ft4^3lb-1/3.  According  to  Szymczak  &  Solomon  (1996), 
Swisdak  published  a  higher  value  for  J  following  reappraisal  of  the  data  presented  by 
Swift  &  Decius  (1947). 

The  period  data  are  central  to  the  calculation  here  and  since  error  bounds  on  the 
estimates  of  7,  e  and  k  are  required,  an  error  estimate  of  the  period  coefficient  K  is 
necessary.  The  data  presented  by  Swift  &  Decius  (1947)  give  K  =  4.36±0.014sft5/6lb-1/3, 
where  the  error  is  the  standard  deviation  of  the  mean.  Converting  to  SI  units  and  retaining 
three  significant  figures  gives  K  =  2.11  ±  O.Olsm^kg"1/3. 

Data  for  the  peak  pressure  experienced  in  the  fluid  surrounding  the  bubble  at  about 
the  time  of  rebound  are  central  to  this  study,  and  the  technique  used  by  Arons  (1948). 
Arons  et  al.  (1948)  presented  such  data  for  the  pressure  measured  at  a  distance  Rp  from 
the  charge  that  satisfied 


Rp  =  w,1^3 /0.352  ft,  Rp  =  w1/3/0.887m,  (4) 

where  is  the  charge  mass  in  pounds  and  w  is  the  charge  mass  is  in  kilograms.  The 
peak  pressure  was  measured  for  charge  masses  0.229  kg  (0.505  lb)  and  1.137  kg  (2.507  lb) 
at  depths  of  76.2  m  (250  ft)  and  152.4  m  (500  ft),  and  charges  of  mass  5.448  kg  (12.01  lb) 
at  a  depth  of  152.4  m  (500  ft). 

Arons  et  al.  (1948)  used  cylindrical  charges  and  measured  the  pressure  at  the  distance 
Rp  at  positions  perpendicular  to  the  cylindrical  surface,  above,  below  and  to  the  side 
of  the  charge.  For  the  largest  charge  of  5.448  kg,  a  peak  pressure  over  hydrostatic  of 
(7.860  ±  0.069)  x  106  Pa  (1140  ±  10  lb/in2)  was  recorded,  with  the  error  the  standard 
deviation  of  the  mean.  Measurements  were  also  made  at  a  distance  Rp  from  the  charge 
along  its  axis.  The  complicating  factor  here  was  that  the  charge  was  detonated  by  a 
blasting  cap  at  one  end  of  the  cylinder.  For  the  5.448  kg  charge  at  a  depth  of  152.4  m, 
the  peak  over-pressures  measured  in  this  configuration  were  8.274  x  106  Pa  (1200  lb/in2) 
at  the  end  closest  the  blasting  cap  and  7.722  x  106  Pa  (1120  lb/in2)  at  the  other  end. 

In  the  methodology  described  in  this  report  it  is  intended  to  provide  error  bounds 
on  the  estimates  of  the  parameters  7,  e  and  k.  A  conservative  approach  is  followed 
and  the  peak  over-pressure  at  the  first  minimum  for  the  5.448  kg  charge  is  taken  to  be 
(7.929  it  0.345)  x  106  Pa  (1150  ±50  lb/in2).  This  range  encompasses  all  the  values  quoted 
above.  Note  also  that  Arons  (1948)  used  a  value  of  8.274  x  106  Pa  (1200  lb/in2)  in  his 
determination  of  the  explosion  bubble  parameters. 

For  the  main  calculation  presented  here,  the  data  for  the  5.448  kg  charge  are  used. 
In  this  case  the  problem  of  accounting  for  the  mass  of  the  detonator  is  minimised.  For 
comparison,  calculations  have  been  performed  using  data  for  the  1.137  kg  charge.  The 
pressure  measurements  for  the  0.229  kg  charge  show  a  much  greater  degree  of  scatter  and 
hence  are  not  used.  Further,  the  objective  of  this  report  is  to  characterise  the  behaviour 
of  larger  quantities  of  explosive. 
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3  Review  of  the  boundary  integral  method 


Boundary  integral  methods  are  well  suited  to  computing  the  motion  of  underwater 
explosion  bubbles,  from  a  time  shortly  after  the  detonation  when  the  shock  wave  has 
travelled  some  distance  from  the  bubble  and  the  fluid  may  be  considered  incompressible. 
The  strategy  adopted  in  studies  of  underwater  explosion  bubble  motion  using  the  boundary 
integral  method  has  been  to  assume  that  the  bubble  is  initially  spherical  in  shape  and  its 
surface  at  rest.  A  high  initial  pressure  within  the  bubble  drives  the  motion.  Although 
these  initial  conditions  are  somewhat  artificial,  the  key  is  to  ensure  that  upon  collapse 
the  bubble  rebounds  at  a  volume  very  close  to  that  at  which  rebound  would  occur  in 
reality.  Thus  the  initial  conditions  are  chosen  to  ensure  that  this  condition  takes  place. 
In  the  very  early  phases  of  the  motion,  just  after  the  shock  wave  is  transmitted  into  the 
water,  buoyancy  and  the  presence  of  nearby  boundaries  affect  the  bubble  motion  little. 
Hence  any  discrepancy  between  the  real  conditions  and  those  used  to  commence  boundary 
integral  calculations  are  not  expected  to  cause  significant  error  in  the  predictions  of  bubble 
behaviour  upon  collapse,  that  being  the  aim  of  the  exercise.  It  is  further  noted  that  Blake 
et  al  (1986)  demonstrated  that  the  growth  phase  of  a  bubble  may  have  a  significant 
effect  on  the  collapse  behaviour  and  jet  formation  process.  Hence  it  is  essential  that  initial 
conditions  are  employed  that  permit  calculation  of  most  of  the  bubble  expansion  phase. 


The  boundary  integral  method  follows  from  the  assumption  that  the  fluid  is  ideal 
and  the  flow  induced  by  the  bubble  motion  irrotational.  A  velocity  potential  <j>  is  then 
introduced  such  that  the  fluid  velocity  is  u  —  V<£.  For  flow  in  a  singly  connected  domain, 
the  potential  is  written  in  the  boundary  integral  form 

c(p)<t>{p)  =  (j^(q)G{p,  q)  -  4>{q)^{p,  9))  dS(q),  (5) 


with 


*>-{£ 


p  e  dft, 
p  e  ft\3ft. 


(6) 


In  these  expressions  p  is  the  position  vector  of  some  point  in  the  flow  field,  G  is  the  Green’s 
function,  ft  the  flow  domain  and  <9ft  its  boundary.  The  normal  to  dCl  is  n,  directed  exterior 
to  ft,  and  d/dn  is  the  normal  derivative.  In  problems  of  bubble  dynamics,  <j)  is  typically 
known  on  free  surfaces,  including  the  bubble  surface,  d<j>/dn  is  known  on  rigid  boundaries 
and  a  solution  of  (5)  is  sought.  Of  particular  relevance  to  the  bubble  is  solution  for  the 
normal  fluid  velocity  d<f>/dn  at  its  boundary. 


The  Bernoulli  equation  is 
dd)  I 

-jfi  +  2  IWP  +  Prv/Pw  +  9(z  ~0=  Poo/pw,  (7) 

where  pw  and  pw  are  the  pressure  and  density  of  the  water  and  g  is  the  gravitational 
acceleration.  The  hydrostatic  pressure  at  the  depth  £  at  which  detonation  occurs  is 
given  by 

Poo  =  Pa  +  Pw9t  =  Pw9(Zo  +  £)>  (8) 

where  the  atmospheric  pressure  pa  =  Pw9i o-  The  Bernoulli  equation  allows  d(f>/dt  to  be 
eliminated  from  the  equation 


d(j>/dt  =  d<j>/dt  +  |V</>|2, 


(9) 
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which  gives  the  rate  of  change  of  </>  following  a  fluid  element.  Suppose  the  pressure  within 
the  bubble  is  specified  independently  of  the  details  of  the  flow  field.  Then  if  (7)  is  evaluated 
at  the  bubble  surface,  pw  is  known  and  (9)  may  be  integrated  in  time  to  determine  the 
evolution  of  <f>  at  the  bubble  surface.  Simultaneous  integration  of  the  equation 

dXjdt  —  u ,  (10) 

describing  the  evolution  of  free  surface  particles  in  the  Lagrangian  sense,  allows  the  motion 
of  the  surface  to  be  computed  in  time.  In  (10)  X  is  the  position  vector  of  a  point  on  the 
bubble  surface  and  u  is  the  fluid  velocity  there.  The  Cartesian  components  of  the  fluid 
velocity  at  the  bubble  surface  are  obtained  from  the  normal  fluid  velocity  component, 
d<fi/dn ,  and  the  tangential  component,  obtained  from  knowledge  of  <j>  over  the  bubble.  In 
accordance  with  (1),  the  pressure  p  within  the  bubble  is 

P  =  Po(V0/V)\  (11) 

where  V  is  the  volume,  Vo  is  its  initial  value  and  po  is  the  initial  value  of  the  pressure  that 
drives  the  bubble  motion  from  rest. 

To  bring  out  the  key  parameters  that  characterise  the  boundary  integral  approach, 
as  introduced  by  Blake  et  al.  (1986)  and  later  Best  &  Kucera  (1992),  the  formulation  is 
transformed  to  dimensionless  variables.  The  distance  scale  is  Rm,  the  maximum  radius  to 
which  the  bubble  would  expand  if  it  remained  spherical,  the  time  scale  is  Rm(pw/Poo )^2 
and  the  pressure  scale  ispoo.  Denoting  dimensionless  quantities  by ',  the  Bernoulli  equation 
becomes 

%  +  ||V^|2  +  aty/VSF  +  -  O  =  1-  (12) 

The  non-dimensional  parameters  that  characterise  the  motion  are  the  ratio  of  specific 
heats  7,  the  buoyancy  parameter 

&  =  (PwdRm/Poo)*^  i  (13) 

and  the  strength  parameter 

a  =  Po/Pco-  (14) 

As  a  buoyant  bubble  collapses  in  an  infinite  fluid  its  added  mass  falls  and  it  accelerates 
upwards.  This  acceleration  causes  a  local  pressure  maximum  to  develop  behind  the  bubble 
and  this  drives  the  motion  of  the  bubble  surface  here  and  causes  formation  of  a  liquid 
jet  that  travels  through  the  bubble.  At  a  given  value  of  the  strength  parameter,  for 
small  buoyancy  parameters  only  part  of  the  jet  may  be  formed  when  the  bubble  achieves 
minimum  volume.  As  the  bubble  re-expands  the  jet  usually  continues  to  travel  through 
the  bubble  (Best  &  Kucera  1992).  For  larger  values  of  the  buoyancy  parameter  the  jet 
may  completely  penetrate  the  bubble  prior  to  it  achieving  minimum  volume.  For  a  given 
buoyancy  parameter  similar  behaviour  is  observed  as  the  strength  parameter  is  varied 
(Best  1991).  For  small  strength  parameters  the  jet  does  not  completely  penetrate  the 
bubble  prior  to  minimum  volume.  As  it  is  increased  jet  formation  occurs  earlier  and 
the  jet  ultimately  penetrates  the  bubble  completely  before  minimum  volume  is  reached. 
When  the  jet  impacts  upon  the  far  side  of  the  bubble  it  evolves  into  toroidal  form,  with 
an  attached  vortex  ring.  The  boundary  integral  method  has  been  extended  by  Best  (1993, 
1994)  to  describe  the  motion  of  toroidal  bubbles. 
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To  proceed  beyond  impact,  Best  connects  the  jet  tip  to  the  other  side  of  the  bubble 
over  a  small  surface  which  forms  a  cut  across  which  the  potential  is  discontinuous  by  A <j>. 
Best  (1993)  has  shown  that  in  an  axisymmetric  geometry 

A 4  =  4>+-  4>~,  (15) 

where  (f)+  and  <£_  respectively  refer  to  the  potential  evaluated  at  the  upper  and  lower 
impact  surfaces  just  before  jet  impact.  Insertion  of  the  cut  allows  the  flow  domain  to  be 
considered  singly  connected  and  (5)  used,  with  special  account  taken  of  the  cut  geometry. 
In  such  a  geometry  the  boundary  integral  expression  for  the  potential  becomes 

c{p)4>{p)  =  (j^(q)G(p,q)  -  <t>(q)^(p,q)j  dS{q)  -  A <p  ^ +(p,q)dS{q ),  (16) 

where  T  denotes  the  cut  surface  and  the  subscript  -f  denotes  evaluation  on  its  upper 
surface.  In  Best’s  original  algorithm  the  surface  T  was  advected  with  the  fluid.  In  a  later 
refinement  (Best  1994)  the  surface  T  is  continuously  deformed  so  that  its  geometry  is 
particularly  simple.  Recently,  toroidal  bubble  shapes  computed  using  Best’s  methodology 
have  been  favourably  compared  with  experimental  data  (Blake,  Tomita  k  Tong  1998, 
Tong,  Schiffers,  Shaw,  Blake  k  Emmony  1999). 

The  pressure  may  be  calculated  anywhere  in  the  flow  field  using  the  Bernoulli  equation 
(7).  The  potential  at  and  around  a  point  is  calculated  using  (5)  or  (16)  and  these  values 
enable  accurate  computation  of  the  time  and  space  derivatives  appearing  in  (7)  using 
simple  finite  difference  approximations. 

Determination  of  the  boundary  integral  parameters  a  and  <5  necessitates  knowledge 
of  the  initial  pressure  po  and  maximum  bubble  radius  Rm-  The  initial  pressure  may  be 
determined  from  assumption  of  the  form  (1)  for  the  variation  of  the  bubble  pressure  and 
the  assumption  that  the  bubble  is  initially  spherical  and  at  rest.  Further,  the  distance  scale 
Rm  is  the  maximum  radius  to  which  the  bubble  would  expand  in  the  absence  of  perturbing 
effects,  that  is  in  the  event  that  the  bubble  remained  spherical.  Thus  to  determine  a  and 
S  in  terms  of  the  physical  parameters  7,  e  and  «,  the  motion  of  a  spherical  bubble  must 
be  considered. 


4  The  spherical  bubble  model 

Suppose  that  a  bubble  undergoes  purely  radial  oscillations  at  depth  £.  The  radius  of 
the  bubble  is  given  by  the  time  dependent  function  R.  Supposing  that  the  origin  is  at  the 
centre  of  the  bubble  the  velocity  potential  for  radial  oscillation  is 

4>=  -R2R/\r\,  (17) 

where  the  dot  denotes  differentiation  with  respect  to  time  and  r  is  the  position  vector  of 
some  point  in  the  flow  field. 

The  usual  way  to  proceed  to  an  equation  of  motion  is  via  energy  conservation.  The 
fluid  kinetic  energy,  T,  is  obtained  as 

T=lpwfv  1 V^'2  dV  =  (18) 
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where  V  denotes  the  volume  occupied  by  the  fluid.  The  potential  energy  of  the  bubble  is 
equal  to  the  work  done  in  expanding  to  radius  R  from  zero  radius  against  the  hydrostatic 
pressure  and  is  equal  to 


(19) 


The  internal  energy  of  the  bubble  contents  is  equal  to  the  work  done  in  compressing  this 
gas  from  zero  pressure,  or  infinite  volume,  to  the  bubble  volume  adiabatically.  Hence  it  is 
equal  to 


(20) 


where  V  =  47ri?3/3  is  the  bubble  volume.  If  w  denotes  the  mass  of  explosive,  considered 
equal  to  the  mass  of  the  detonation  gas,  then  from  (1)  the  pressure  along  an  adiabatic 
characterised  by  k  is 

p  =  k(w/VF,  (21) 


so 


\3  /  7  —  1 


(22) 


Inr  consideration  of  (21)  it  is  pertinent  to  note  that  photographic  records  of  underwater 
explosions  suggest  that  not  all  the  explosive  mass  is  converted  to  gas  upon  detonation.  The 
photographs  show  what  appear  to  be  solid  remnants  of  the  detonation,  which  sometimes 
obscure  the  image  of  the  bubble  surface.  It  may  therefore  be  appropriate  to  assume  that 
a  fixed  proportion  of  the  mass  is  converted  to  gas.  This  would  cause  the  appearance  of  an 
additional  constant  within  (21).  However,  this  constant  may  simply  be  incorporated  into 
a  redefined  value  of  k  and  the  analysis  that  follows  is  unchanged.  Given  the  assumptions 
made  in  modelling  the  water  and  explosion  contents  it  would  be  inconsistent  to  assume 
that  anything  but  a  constant  proportion  of  the  explosive  mass  were  converted  to  gas  upon 
detonation. 


The  equation  of  motion  is  obtained  by  equating  the  sum  of  these  energies  to  the  total 
energy  in  the  bubble  oscillation,  which  is  equal  to  ew.  To  simplify  proceedings,  appropriate 
time  and  length  scales  are  introduced.  The  distance  scale  is 


Rsc  —  ^ 


3  ew 


(23) 


and  the  time  scale  is 

r— <«> 

Note  that  these  are  different  from  the  scales  introduced  in  the  previous  section.  The 
equation  for  the  non-dimensional  bubble  radius  r  is  then 


r3r2  +  rs  +  fir  3^7  ^  =  1, 


(25) 


where 


M  = 


(7  —  i)*7 


(26) 


At  the  maximum  and  minimum  radii,  r  =  0,  so  they  are  given  as  the  solution  of 


r3  4*  fir  3^7  ^  —  1  =  0. 


(27) 
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Denoting  them  as  rmax  and  rm\n  they  are  given  approximately  when  p  C  1  as 


rmin  ~  P 


l/(3(7-l)) 

5 


(28) 


and 

W«l-/i/3+(2-37)M2/9-  (29) 

Given  these  approximate  values  the  solutions  may  be  found  numerically  using  bisection, 
say. 

The  solutions  of  (25)  have  the  characteristic  that  they  are  symmetric  about  the  maxi¬ 
mum  radius.  Hence  the  non-dimensional  period  of  oscillation  is  given  by 


=  2/ 

“  T min 


-3/2 


(1  —  —  ^r“3(7-*i))i/2 


dr. 


(30) 


Evaluation  of  (30)  is  performed  numerically,  using  Gaussian  quadrature,  with  the  singu¬ 
larities  in  the  integrand  at  both  ends  of  the  integration  range  subtracted  and  integrated 
analytically.  In  the  limiting  case  /x  =  0,  the  integral  may  be  evaluated  analytically  in 
terms  of  the  Gamma  function  and  yields  r  =  v/3/7rF(5/6)r(2/3)  =  1.4937. 

For  future  reference  it  is  worthwhile  to  write  down  expressions  for  the  maximum  and 
minimum  radii  and  the  oscillation  period.  The  period  is 


/  3  \ 5/6 /  e  \V»  tpV* 

V  2ff  /  V  2npw  )  (4  +  $o)5/6  T’ 

the  maximum  radius  is 


and  the  minimum  radius  is 


(31) 

(32) 

(33) 


Using  the  minimum  radius,  po  may  be  evaluated  from  (1)  and  the  strength  parameter 
a  is  given  in  terms  of  the  explosive  charge  parameters  by. 

<X  =  «P3o_1/(er min)7  =  (7  -  1  Warnin'  C34) 

The  buoyancy  parameter  is 

*  =  (PwgRm/Poo)1'2  =  (Rm/(i  +  6>))1/2.  (35) 


For  later  comparison  an  expression  is  required  for  the  pressure  experienced  in  the 
fluid  at  bubble  minimum.  In  dimensional  variables,  the  pressure  is  given  by  the  Bernoulli 
equation  (7)  as 


Pw/Pw  =  Poo/Pw  -  ^  -  ||W|2, 


(36) 
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noting  that  the  term  involving  the  gravitational  acceleration  is  not  consistent  with  the 
spherical  bubble  assumption.  From  (17)  it  follows  that 

^  =  -(2RR?  +  R2R)/\r\,  (37) 

and 

|V<£|2  =  R4R2/\r\4.  (38) 

The  maximum  pressure  occurs  at  bubble  minimum  when  R  =  0  so 

Pma»/Poo  -  1  =  (39) 

3  \r\ 

where  the  distance  and  time  scales  (23)  and  (24)  have  been  used  and  rmjn  is  the  second 
time  derivative  of  the  non-dimensional  radius  at  minimum  volume.  Differentiating  (25) 
with  respect  to  time,  dividing  through  by  r  and  then  setting  it  equal  to  zero  yields 

2  =  rm.n((7  _  1)/ir-37  _  1})  (40) 

where  (27)  has  been  used.  Noting  the  definition  of  the  strength  parameter  a  in  (34),  the 
peak  pressure  at  bubble  minimum  may  be  written  as 


PmaxM/Poo  -  1  =  (a  -  l)rmin/|r|.  (41) 

Note  that  in  this  expression  |r|  is  non-dimensional  with  respect  to  i^c* 


5  Determination  of  the  parameters 

The  first  point  to  be  made  with  respect  to  a  determination  of  the  three  bubble  pa¬ 
rameters  7,  e  and  k  following  the  technique  of  Arons  (1948)  is  that  only  two  pieces  of 
information,  namely  equation  (2)  and  peak  pressure  data  are  available.  Thus  the  strat¬ 
egy  to  be  followed  is  to  determine  the  three  parameters  so  that  the  theoretical  bubble 
period  differs,  in  a  least  squares  sense,  from  (2)  by  as  little  as  possible,  and  that  the  peak 
pressure  at  bubble  minimum  is  in  agreement  with  the  experimental  data.  The  essential 
difference  between  this  work  and  that  of  Arons  is  that  here  the  restrictive  spherical  bub¬ 
ble  assumption  is  removed  through  use  of  the  boundary  integral  method  to  compute  the 
bubble  motion  and  peak  pressure. 

That  the  bubble  period  data,  as  opposed  to  maximum  radius  data,  should  be  central  to 
any  determination  of  the  bubble  parameters  stems  from  the  fact  that  this  characteristic  of 
the  explosion  event  may  be  easily  measured  from  a  pressure  record.  Data  for  the  maximum 
bubble  radius  requires  that  high  speed  photography  be  employed  in  such  a  manner  that 
reliable  distance  measurements  may  be  undertaken.  In  addition,  it  is  often  found  in  high 
speed  photography  of  underwater  explosion  events  that  the  bubble  surface  is  not  smooth, 
a  feature  that  introduces  error  into  determinations  of  the  maximum  radius. 

Before  embarking  on  the  computations  it  is  noted  that  deformation  of  the  bubble  from 
spherical  shape  does  not  affect  the  bubble  period  by  an  amount  that  is  significant  in  the 
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context  of  these  calculations.  A  posteriori  calculations  show  that  for  a  charge  of  mass 
5.448  kg  over  the  range  of  depths  from  10  m  -  150  m,  the  period  as  computed  using  the 
boundary  integral  method  differs  from  the  spherical  bubble  value  by  at  most  0.67%.  For 
depths  of  50  m  and  below  the  deviation  is  at  most  0.15%.  This  deviation  is  negligible 
given  uncertainty  in  the  experimental  data.  Hence  the  spherical  bubble  expression  for  the 
period  is  used. 

To  proceed,  compare  the  expressions  of  (2)  and  (31).  It  is  assumed  that  the  three 
bubble  parameters  are  constants  and  not  dependent  upon  the  charge  mass  or  depth  of 
detonation.  In  this  event,  equating  (2)  and  (31)  indicates  that  the  non-dimensional  bubble 
period  r  is  constant.  It  is  clear  from  (30)  that  r  is  dependent  upon  depth  via  its  dependence 
upon  p.  Consider  figure  1,  which  displays  the  non-dimensional  bubble  period  as  a  function 
of  p  for  different  values  of  7.  This  figure  has  been  generated  following  Swift  &  Decius 
(1947)  and  Arons  (1948).  Values  of  p  that  are  realised  in  practice  generally  range  over 
the  interval  (0,0.4),  a  fact  that  is  confirmed  a  posteriori ,  once  the  parameters  have  been 
estimated.  Arons  used  such  a  figure  to  select  a  value  of  7  such  that  the  variation  of  the 
non-dimensional  period  would  be  minimal,  following  comparison  of  (2)  and  (31).  Arons 
chose  a  value  of  1.25.  The  value  of  5/4  also  renders  the  exponents  in  (25)  and  (26)  easily 
simplified. 


Figure  1:  The  non-dimensional  bubble  period  r  as  a  function  of  p  for  a  selection  of  values 

0/7* 

This  technique  is  to  be  followed  here,  but  with  greater  precision.  It  is  supposed  that  (2) 
is  valid  over  a  range  of  depths,  here  taken  to  be  0  - 150  m,  in  accordance  with  comments  by 
Arons  (1948).  Given  values  for  the  parameters  7,  e  and  k,  this  range  can  be  equivalently 
considered  a  range  of  p.  Suppose  it  takes  on  values  in  the  range  (pmm>  /Wx)  for  some 
maximum  and  minimum  values.  Denote  A p  =  pm*x  —  Pmin-  Over  this  range  the  mean 


10 


DSTO-RR-0238 


non-dimensional  period  is 

max 

f  (7)  =  A/x  1  /  r(/x,7)d/A,  (42) 

**  Mmin 

and  the  mean  square  deviation  is 

-  1  /“Mmax 

At2 (7)  =  Ap  1  /  (r(p,7)  -f(j))2dn, 

*  min 
_  /'Aimax 

=  Ap-1  /  r2(/x,  7)  -  f 2(7).  (43) 

•'Mmin 


The  procedure  adopted  is  to  find  the  value  of  7  such  that  Ar2  is  a  minimum.  It  should 
be  noted  that  7  is  restricted  to  lie  in  the  regime  of  physical  values  for  gases,  namely 
(1, 5/3].  In  the  absence  of  values  for  7,  e  and  k  the  initial  choice  for  the  range  (/Umin,  /zmax) 
is  (0.0,0.35).  The  integrals  of  (42)  and  (43)  are  evaluated  using  Gaussian  quadrature  and 
the  minimum  is  determined  using  a  golden  section  search  (Press  et  al  1992). 

Denote  by  r  the  mean  value  of  r  at  the  value  of  7  determined  by  minimisation  of  (43). 
The  bubble  energy  per  mass  of  explosive  is  then  given  by  equating  (2)  and  (31)  and  is 

e  =  2np(2g/3)5/2(K/f)3.  (44) 

The  remaining  part  of  the  calculation  involves  a  determination  of  ac  in  order  to  reconcile 
theoretical  predictions  of  the  peak  pressure  upon  bubble  rebound  with  the  experimental 
data  of  Arons  et  al  (1948).  Following  estimation  of  ac,  a  more  reliable  estimate  of  the  range 
of  values  of  \i  may  be  made.  The  procedure  for  estimating  the  three  bubble  parameters 
is  then  repeated  to  yield  refined  estimates.  Given  the  accuracy  of  the  experimental  data 
and  numerical  limitations,  the  procedure  is  iterated  only  once. 

It  is  appropriate  to  comment  that  two  parameters  are  essentially  determined  on  the 
basis  of  the  bubble  period  data.  It  must  be  emphasised  that  this  characteristic  of  the 
bubble  event  is  by  far  the  easiest  to  measure  and  it  is  therefore  appropriate  that  it  be 
utilised  to  maximum  extent  in  determination  of  the  physical  parameters  that  characterise 
the  bubble  contents. 

To  proceed  to  a  determination  of  ac  requires  that  pmax(Rp)  he  determined  as  a  function 
of  ac,  following  estimation  of  7  and  e.  Given  that  a  computation  of  the  whole  bubble  motion 
using  the  boundary  integral  method  is  required  to  determine  this  pressure  for  each  value  of 
ac,  the  technique  used  was  designed  to  minimise  the  number  of  computations.  The  required 
pressure  was  evaluated  for  ac  equal  to  1.0  x  105, 1.125 x  105,  1.25  x  105,  1.375  x  105,  1.5  x  105 
and  1.625  x  105,  values  that  were  found  to  well  cover  the  region  where  the  solution  was 
found.  A  cubic  spline  was  fitted  to  the  data  to  interpolate  a  smooth  curve  and  a  bisection 
method  used  to  find  the  value  of  k  at  which  the  theoretical  prediction  of  pressure  coincided 
with  the  experimental  data. 

In  order  to  provide  some  estimate  of  the  accuracy  of  the  parameters  determined,  the 
procedure  was  repeated  using  values  for  the  coefficient  K  and  the  peak  pressure  that  are  at 
the  extremes  of  the  ranges  defined  by  the  error  estimates  quoted  in  section  2.  Calculating 
over  all  combinations  of  extreme  values  provides  extreme  estimates  of  7,  e  and  ac  and  hence 
estimates  of  error. 
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To  perforin  the  calculations  presented  here  requires  values  to  be  fixed  for  the  fluid 
density,  atmospheric  head  and  gravitational  acceleration.  These  values  are  taken  to  be  as 
follows: 

pw  =  1027.6  kg  m-3, 

&  =  10.05  m, 
g  —  9.81ms-2. 

This  value  for  the  fluid  density  is  that  of  sea  water  at  5°  C.  The  boundary  integral  method 
computes  pressures  scaled  with  respect  to  the  hydrostatic  pressure  at  the  depth  of  deto¬ 
nation.  With  this  scaling,  the  peak  pressure  recorded  by  Arons  et  al  (1948)  for  a  charge 
of  mass  5.448  kg  at  a  depth  of  152.4  m,  quoted  in  section  2,  is  equal  to  5.842  ±  0.210. 

The  values  obtained  for  the  parameters  characterising  the  bubble  contents  are: 

7  =  1-23, 

6  =  (2.08  ±  0.03)  x  106J  kg-1, 

k  =  (1.46  ±  0.05)  x  105, 

noting  that  to  three  significant  figures  7  is  errorless. 

For  comparison,  the  calculation  was  also  performed  using  peak  pressure  data  for  the 
1.137  kg  charge  detonated  at  a  depth  of  76.2  m.  The  peak  over-pressure  recorded  by 
Arons  et  al  (1948)  was  7.033  x  106  Pa  (1020  lb/in2).  The  values  obtained  are:  7  =  1.23, 
e  =  2.09  x  106  Jkg“\  k  =  1.53  x  105.  Only  the  value  for  k  is  outside  the  quoted  error 
range,  indicating  the  sensitivity  of  this  parameter  to  the  pressure  data. 

The  values  determined  above  are  close  to  those  given  by  Arons  (1948): 

7  =  1.25,  €  =  2.05  x  106  Jkg“\  k  =  1.45  x  105. 

That  the  values  for  7  and  e  should  be  close  follows  from  their  determination  primarily 
using  period  data,  in  conjunction  with  the  observation  that  deformation  of  the  bubble  from 
spherical  shape  during  collapse  has  negligible  effect  upon  the  period.  Closer  investigation 
of  the  effect  of  bubble  deformation  upon  the  pressure  field  surrounding  the  bubble  shows 
the  close  agreement  of  k  to  be  somewhat  fortuitous. 


6  Variation  of  the  pressure 

Having  determined  the  parameters  characterising  the  bubble  contents,  it  is  interesting 
to  compute  the  pressure  field  around  the  bubble,  and  particularly  the  maximum  value  as 
the  bubble  rebounds  for  the  first  time.  The  calculations  presented  here  are  of  the  pressure 
at  a  distance  from  the  charge  that  satisfies  (4)  and  these  are  compared  with  the  spherical 
bubble  expression  (41). 

First,  the  pressure  has  been  calculated  at  a  radius  satisfying  (4),  for  a  range  of  charge 
masses  between  1  kg  and  1000  kg,  at  a  depth  of  152.4  m,  which  replicates  the  conditions 
under  which  the  experimental  data  used  here  were  compiled.  The  results  are  shown  in 
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figure  2.  The  unevenness  of  this  curve  arises  from  effects  associated  with  numerical  differ¬ 
entiation  used  in  computation  of  the  pressure.  Also  shown  is  the  expression  obtained  from 
the  spherical  bubble  model.  In  this  model,  for  a  given  depth,  the  pressure  is  constant  for 
w  and  Rp  satisfying  (4).  The  results  obtained  using  the  boundary  integral  method  show 
significantly  more  complex  behaviour  of  the  peak  pressure  and  this  is  due  to  translation 
and  deformation  of  the  bubble  from  spherical  shape. 


Figure  2:  The  non-dimensional  peak  pressure  experienced  at  Rp  due  to  bubble  rebound, 
as  a  function  of  w1/3.  The  depth  is  154.2  m.  The  smooth  curve  is  the  prediction  of  the 
spherical  bubble  model 

For  small  charge  masses  the  peak  pressure  is  close  to  the  value  predicted  by  the  spher¬ 
ical  model.  As  the  charge  mass  increases  the  peak  pressure  falls  and  reaches  a  local 
minimum,  after  which  it  begins  to  increase  rapidly.  Inspection  of  the  bubble  shapes  dur¬ 
ing  the  collapse  and  rebound  demonstrates  features  that  correlate  with  these  features  of 
the  peak  pressure  and  hence  assist  in  interpreting  the  behaviour.  Consider  the  case  where 
w1/3  =  1.76,  which  corresponds  to  the  data  presented  by  Arons  et  al .  (1948)  and  used 
here.  The  bubble  collapse  and  rebound  are  shown  in  figure  3,  the  bubble  being  spherical 
to  a  high  degree  of  approximation  during  the  first  expansion,  which  is  not  shown.  The 
critical  observation  is  that  the  bubble  is  close  to  spherical  as  it  rebounds.  A  flattening  of 
the  lower  surface  is  evident  at  minimum  volume  and  is  a  prelude  to  jet  formation  upon 
bubble  rebound.  As  the  bubble  rebounds  the  jet  develops  and  travels  through  the  bubble, 
as  first  documented  by  Best  &  Kucera  (1992).  Once  the  jet  completely  penetrates  the 
bubble,  an  event  referred  to  as  jet  impact,  it  evolves  into  toroidal  form,  as  shown  in  figure 
4.  There  is  a  weak  circulation  around  the  torus.  As  the  bubble  continues  to  re-expand  the 
thread  of  fluid  flowing  through  the  torus  thins  and  it  appears  that  the  bubble  will  resume 
a  singly  connected  topology.  Computation  of  the  subsequent  behaviour  has  been  given  by 
Best  (1994). 

Given  the  degree  to  which  the  bubble  may  be  approximated  as  spherical  at  the  time  of 
rebound,  it  is  no  surprise  that  the  boundary  integral  method  and  spherical  bubble  model 
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axe  in  approximate  agreement  as  to  the  peak  pressure.  As  the  charge  mass  increases  so 
does  the  degree  of  jet  formation  at  bubble  minimum.  The  quantity  of  fluid  kinetic  energy 
at  the  minimum  therefore  increases  and  as  a  consequence  of  energy  conservation  the  bubble 
does  not  collapse  to  as  small  a  volume  as  it  would  if  it  remained  spherical  and  stationary. 
This  reasoning  is  consistent  with  the  evidence  of  figure  5  which  shows  the  variation  with 
charge  mass  of  the  bubble  volume  at  the  first  minimum.  The  time  derivative  of  the  velocity 
potential  around  the  time  of  first  minimum  is  dominated  by  a  term  that  reflects  the  rate 
of  change  of  bubble  volume.  As  the  minimum  volume  increases  with  charge  mass,  its  time 
derivative  and  that  of  the  velocity  potential  decreases.  The  non-dimensional  peak  pressure 
at  Rp  therefore  decreases  with  charge  mass. 
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Figure  3:  The  collapse  (a)  and  rebound  (b)  of  the  bubble  created  by  detonating  5.44 8  k 9 
TNT  at  a  depth  of  154.2  m.  The  distance  scale  is  non-dimensional  with  respect  to  the 
maximum  bubble  radius ,  equal  to  1.04  The  times,  in  milliseconds,  corresponding  to  the 
bubble  shapes  are:  (a)  34.92,  48.01,  52.06,  52.95,  53.33,  (b)  53.58,  53.82,  54.11 ,  H-tf, 
54.16,  54.98. 


-1  00  -0.25  0.50  -1.00  -0.25  0.50  -1.00  -0.25  0.50 

r  r  r 


Figure  4:  The  motion  of  the  toroidal  bubble  formed  during  the  rebound  of  the  bubble  created 
by  detonating  5.448  kg  TNT  at  a  depth  of  154.2  m.  The  distance  scale  is  non-dimensional 
with  respect  to  the  maximum  bubble  radius,  equal  to  1.04  m.  The  times,  in  milliseconds, 
corresponding  to  the  bubble  shapes  are:  (a)  54.98,  (b)  55.70,  (c)  56.15. 

As  the  charge  mass  is  increased  further  the  stage  is  reached  where  the  bubble  evolves 
into  toroidal  form  before  minimum  volume  is  reached.  An  example  is  illustrated  in  figures 
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6  and  7,  for  the  case  =  7.  Under  such  circumstances  the  bubble  at  minimum  volume 
has  an  attached  vortex  ring  and  the  peak  pressure  in  this  regime  is  an  increasing  function 
of  the  charge  mass.  Notice  that  the  local  minimum  in  the  peak  pressure  curve  of  figure 
2  corresponds  closely  with  that  charge  mass  at  which  the  bubble  becomes  toroidal  at 
minimum  volume.  It  appears  that  this  extremum  in  the  peak  pressure  is  a  signature  of 
this  change  in  physical  behaviour  from  the  bubble  being  singly  connected  at  minimum 
volume  to  being  a  toroidal  bubble  with  an  attached  vortex  ring  at  this  time. 


Figure  5:  The  volume  of  the  bubble  at  the  first  minimum ,  as  a  function  of  The 

depth  of  detonation  is  154-2  m.  The  volume  is  non-dimensional,  with  the  distance  scale 
the  maximum  bubble  radius. 

In  this  context  it  is  worthwhile  recalling  computations  of  the  pressure  field  in  the  fluid 
before  and  after  toroidal  bubble  formation  as  presented  by  Best  (1993).  These  demonstrate 
that  before  jet  impact  there  is  a  local  maximum  in  the  pressure  field  located  behind  the 
bubble  and  driving  the  jet  through  it.  After  jet  impact  the  location  of  this  local  maximum 
changes  to  a  position  ahead  of  the  bubble  and  the  value  of  the  maximum  is  enhanced. 
The  action  of  the  peak  pressure  here  is  to  decelerate  the  rush  of  fluid  through  the  torus 
and  drive  the  circulatory  motion  of  fluid  around  the  torus.  The  local  minimum  in  the 
peak  pressure  as  a  function  of  charge  mass  noted  in  figure  2  clearly  reflects  the  change  in 
character  of  the  pressure  field  brought  about  by  toroidal  bubble  formation.  For  smaller 
charge  masses  the  peak  pressure  at  rebound  reflects  the  pressure  field  associated  with  a 
singly  connected  bubble  that  exhibits  a  local  maximum  behind  the  bubble.  For  larger 
charge  masses  the  peak  pressure  at  bubble  rebound  reflects  the  pressure  field  associated 
with  a  toroidal  bubble  that  exhibits  a  local  maximum  ahead  of  the  bubble. 

It  is  instructive  to  interpret  the  changing  behaviour  of  the  bubble  with  charge  mass 
in  terms  of  the  non-dimensional  parameters  introduced  in  conjunction  with  the  boundary 
integral  method.  The  strength  parameter  a  is  independent  of  charge  mass,  so  plays  no 
role  in  interpretation  of  the  present  set  of  results.  The  key  parameter  is  the  buoyancy 
parameter,  given  by  (35).  At  a  given  depth  5  varies  as  w 1/6.  The  variation  of  bubble  be¬ 
haviour  with  buoyancy  parameter  has  been  systematically  investigated  by  Best  &  Kucera 
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Figure  6:  The  collapse  of  the  bubble  created  by  detonating  343  kg  TNT  at  a  depth  of  154-2 
m.  The  distance  scale  is  non-dimensional  mth  respect  to  the  maximum  bubble  radius ,  equal 
to  4.13  m.  The  times ,  in  milliseconds,  corresponding  to  the  bubble  shapes  are:  136.6,  188.0, 
204.7,  208.5,  210.1,  211.1,  211.9,  213.1. 

(1992).  The  key  role  that  buoyancy  plays  in  the  jet  formation  process  is  to  cause  the  rapid 
acceleration  upwards  of  the  approximately  spherical  bubble  as  it  begins  to  collapse.  The 
consequence  of  this  is  development  of  a  localised  pressure  maximum  behind  the  bubble 
which  drives  the  jet  through  the  bubble. 

As  the  buoyancy  force,  as  quantified  by  <5,  increases,  jet  formation  occurs  earlier.  For 
small  charge  masses  this  causes  more  energy  to  manifest  itself  in  the  form  of  fluid  kinetic 
energy  at  the  time  of  minimum  volume,  causing  an  increase  in  the  minimum  volume 
and  consequent  downward  trend  in  peak  pressure.  As  the  buoyancy  parameter  increases 
further  a  stage  is  reached  where  jet  penetration  coincides  with  the  bubble  minimum.  As 
6  increases  further  jet  penetration  occurs  earlier  and  rebound  occurs  after  the  bubble  has 
evolved  into  toroidal  form.  In  this  regime,  the  fluid  motion  associated  with  the  circulation 
around  the  torus  causes  the  peak  pressure  to  increase  with  6. 

An  additional  set  of  calculations  has  been  performed  for  the  maximum  pressure  at  Rp 
due  to  bubble  rebound  in  the  case  where  the  depth  is  76.2  m,  which  corresponds  to  the 
other  depth  at  which  Arons  et  al  (1948)  gathered  pressure  data.  The  maximum  pressure 
as  a  function  of  charge  mass  is  shown  in  figure  8.  This  figure  also  displays  the  value  given 
by  the  spherical  bubble  model  as  the  constant  curve.  There  are  two  aspects  to  note  when 
comparing  this  figure  to  figure  2,  which  corresponds  to  a  depth  of  152.4  m.  The  strength 
parameter  a  is  greater  at  the  smaller  depth  and  for  a  given  charge  mass  the  buoyancy 
parameter  varies  as  (f  +  £o)~2/3,  to  leading  order.  Hence  a  smaller  depth  yields  a  larger 
buoyancy  parameter.  This  is  in  accordance  with  the  observation  that  a  given  mass  of 
explosive  will  form  a  larger  bubble  at  smaller  depth  owing  to  the  reduced  hydrostatic 
pressure,  with  the  buoyancy  force  consequently  larger. 
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Figure  7:  The  motion  of  the  toroidal  bubble  formed  following  collapse  of  the  bubble  created 
by  detonating  343  kg  TNT  at  a  depth  of  154-2  m.  The  distance  scale  is  non-dimensional 
with  respect  to  the  maximum  bubble  radius ,  equal  to  4.13  m.  The  times,  in  milliseconds, 
corresponding  to  the  bubble  shapes  are :  (a)  213.1,  (b)  214.2,  (c)  215.9,  (d)  218.6,  (e) 
222.2,  (f)  226.6. 


Figure  8:  The  non-dimensional  peak  pressure  experienced  at  Rp  due  to  bubble  rebound, 
as  a  function  of  wl/3.  The  depth  is  76.2  m.  The  smooth  curve  is  the  prediction  of  the 
spherical  bubble  model. 
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In  this  example,  the  bubble  for  which  w1/3  =  1.5  reached  minimum  volume  just  prior 
to  the  transition  to  a  toroidal  geometry,  and  for  w =  2.0  the  bubble  became  toroidal 
just  prior  to  achieving  minimum  volume.  Hence  the  local  minimum  of  the  peak  pressure 
curve  again  corresponds  with  the  coincidence  of  jet  impact  and  minimum  volume.  That 
this  correspondence  should  occur  at  a  much  smaller  charge  mass  is  due  to  the  two  effects 
noted  above  for  the  strength  parameter  and  buoyancy  parameter.  The  strength  parameter 
provides  some  measure  of  the  rate  of  change  of  bubble  volume,  and  hence  added  mass,  as 
the  bubble  collapses.  The  greater  the  rate  of  change  the  larger  the  upwards  acceleration  of 
the  bubble.  The  increased  buoyancy  parameter  acts  similarly,  causing  earlier  jet  formation. 
The  additional  feature  evident  from  these  data,  compared  with  those  of  figure  2,  is  the 
behaviour  of  the  maximum  pressure  as  the  charge  mass  is  increased.  Once  the  bubble 
becomes  toroidal  before  minimum  volume  is  reached,  the  peak  pressure  at  Rp  increases 
with  w.  However,  a  maximum  is  achieved  and  thereafter  it  falls. 

A  final  series  of  computations  was  performed  to  explore  the  variation  of  peak  pressure 
at  Rp  for  a  given  charge  mass,  as  the  depth  varies.  The  data  axe  shown  in  figure  9  for 
w  =  3.375  kg.  The  pressure  scale  here  is  atmospheric  pressure.  Also  shown  is  the  value 
predicted  by  the  spherical  model,  which  displays  a  monotonic  increase  with  depth.  At 
the  larger  depths  the  two  results  are  in  close  agreement.  This  reflects  the  fact  that  both 
the  strength  and  buoyancy  parameters  decrease  with  depth  and  the  bubble  therefore  is 
approximately  spherical  at  the  time  of  the  first  minimum.  As  the  depth  decreases,  both 
increase  and  cause  jet  formation  to  occur  earlier.  The  coincidence  of  the  transition  from 
simply  connected  to  toroidal  geometry  and  minimum  volume  occurs  at  a  depth  between  70 
m  and  65  m,  and  this  corresponds  to  a  local  minimum  of  the  peak  pressure  variation.  As 
the  depth  decreases  towards  this  extremum,  the  greater  extent  of  jet  formation  at  minimum 
volume  results  in  increased  fluid  kinetic  energy  at  the  expense  of  bubble  internal  energy 
and  hence  the  minimum  volume  increases.  The  peak  pressure  is  thus  below  the  spherical 
bubble  prediction.  As  the  bubble  becomes  toroidal  at  minimum  volume  the  nature  of  the 
flow  field  manifests  itself  as  an  increasing  peak  pressure.  As  the  parameters  a  and  6  that 
quantify  jet  formation  increase  further  with  decreasing  depth  the  peak  pressure  reaches  a 
maximum  and  thereafter  decreases.  This  trend  is  as  noted  in  figure  8. 


To  obtain  further  insight  into  the  factors  affecting  the  variation  of  peak  pressure  with 
depth,  it  is  useful  to  compute  the  peak  pressure  caused  by  bubble  rebound  in  a  spherical 
bubble  model  in  which  bubble  translation  is  taken  into  account.  In  this  case  the  fluid 
velocity  potential  is 

F?R  S?Ue,.  (r-») 

+  ~~  2|r  —  *|3  '  (45) 

where  ez  is  a  unit  vector  in  the  vertical  direction,  z  =  zez  is  the  position  vector  of  the 
bubble  centroid  and  U  =  z  is  the  velocity  of  the  bubble.  Consideration  of  energy  and 
momentum  conservation  yield  the  following  equations  of  motion  (Best  1991): 


RR  +  \k2  -  \u 2  =  ^-(Ro/Rf1  +9*-—,  (46) 

2  4  pw  pw 

£( R3U )  =  2  gR3.  (47) 

at 

In  these  expressions  it  is  assumed  that  at  t  =  0  the  bubble  centroid  is  located  at  z  =  0 
and  the  bubble  is  stationary  so  U  =  0.  These  equations  may  be  solved  numerically  using 
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Figure  9:  The  non-dimensional  peak  pressure  experienced  at  Rp  due  to  bubble  rebound, 
as  a  function  of  detonation  depth.  The  charge  mass  is  3.375  kg  and  the  pressure  scale  is 
atmospheric  pressure.  The  smooth  solid  curve  is  the  prediction  of  the  stationary  spherical 
bubble  model.  The  dashed  curve  is  the  prediction  of  the  translating  spherical  bubble  model. 


standard  techniques  and  use  of  (45)  in  the  Bernoulli  equation  (7)  permits  computation  of 
the  pressure  field  in  the  fluid. 

The  peak  pressure  at  bubble  rebound  as  a  function  of  depth  is  shown  as  the  dashed 
curve  in  figure  9.  Compared  with  the  stationary  spherical  bubble  model,  the  translating 
spherical  bubble  model  shows  a  decrease  in  the  peak  pressure  as  depth  decreases.  This 
trend  is  in  accordance  with  the  variation  computed  using  the  boundary  integral  method. 
At  the  greater  depths,  it  has  been  established  that  the  bubble  remains  spherical  at  the 
first  minimum,  to  a  high  degree  of  approximation.  In  this  regime  the  translating  spherical 
bubble  model  might  be  expected  to  well  describe  the  behaviour. 

It  is  instructive  to  compute  the  maximum  centroid  velocity  during  the  first  bubble  col¬ 
lapse  and  early  rebound  phase  using  both  the  boundary  integral  method  and  the  spherical 
bubble  model.  The  results  are  shown  in  figure  10.  In  this  figure  the  velocity  is  dimension¬ 
less  with  respect  to  that  computed  with  the  boundary  integral  method  at  the  maximum 
depth  of  200  m.  In  this  context,  note  that  the  velocity  scale  varies  as  (£  +  £q)V2,  hence 
the  requirement  to  fix  a  depth  to  scale  all  results. 

These  results  are  quite  interesting  and  assist  in  interpreting  the  behaviour.  They 
show  that,  even  for  quite  large  depths,  when  the  jet  has  only  partially  formed  at  bubble 
minimum,  the  upwards  velocity  of  the  non-spherical  bubble  is  noticeably  less  than  the 
value  predicted  using  the  spherical  model.  In  addition,  figure  9  shows  that  the  pressure 
predicted  using  the  boundary  integral  method  and  translating  spherical  bubble  model 
begins  to  deviate  from  that  predicted  using  the  stationary  spherical  bubble  model  at  a 
depth  of  around  80-100  m.  At  this  depth  the  translating  spherical  bubble  models  predicts 
an  upwards  velocity  about  twice  that  given  by  the  boundary  integral  method. 
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Figure  10:  The  maximum  centroid  velocity  during  the  collapse  phase,  as  a  function  of 
depth ,  computed  using  the  boundary  integral  method.  The  charge  mass  is  3.375  kg  and 
the  velocity  scale  is  the  value  predicted  at  a  depth  of  200  m.  The  smooth  curve  represents 
computations  using  the  translating  spherical  bubble  model. 


Figure  11:  The  motion  of  the  toroidal  bubble  formed  following  collapse  of  the  bubble  created 
by  detonating  3.375  kg  TNT  at  a  depth  of  20  m.  The  distance  scale  is  non-dimensional 
with  respect  to  the  maximum  bubble  radius,  equal  to  1.64  The  times ,  in  milliseconds, 
corresponding  to  the  bubble  shapes  are:  (a)  185.8,  (b)  186.1 ,  (c)  186.3,  (d)  186.7,  (e) 
187.6,  (f)  188.2. 
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If  the  deviation  of  the  pressure  from  that  predicted  using  the  stationary  spherical  bub¬ 
ble  model  was  due  to  translation  alone  then  we  would  expect  the  centroid  velocities  given 
by  the  boundary  integral  method  and  translating  spherical  bubble  model  to  be  similar. 
The  fact  that  they  differ  by  a  factor  of  about  two  demonstrates  that  the  deviation  of  the 
pressure  from  the  stationary  spherical  bubble  model  prediction  is  not  due  to  translation 
alone  but  a  combination  of  translation  and  deformation  from  spherical  shape. 

As  the  depth  decreases  further,  the  translating  spherical  bubble  model  predicts  a  peak 
pressure  that  continues  to  fall,  whereas  the  boundary  integral  method  calculation  shows 
that  the  pressure  has  a  local  minimum  at  the  depth  where  minimum  volume  and  the 
transition  to  a  toroidal  topology  coincide.  This  divergence  in  predictions  indicates  that 
the  variation  of  peak  pressure  with  depth  reflects  the  complex  flow  field  of  the  toroidal 
bubble  with  attached  vortex  ring.  It  also  illustrates  that  in  this  regime  the  stationary 
spherical  bubble  model  is  superior  to  that  which  includes  translation  in  predicting  the 
pressure  in  the  flow  field. 

The  variation  of  peak  pressure  with  depth,  as  computed  using  the  boundary  integral 
method  and  illustrated  in  figure  9,  shows  three  distinct  regimes,  separated  by  a  local 
maxima  and  minima.  It  might  be  suggested  that  these  correspond  to  different  regimes  of 
bubble  behaviour.  The  discussion  above  suggests  that  at  large  depth  values  the  significant 
attribute  of  the  bubble  behaviour  is  that  jet  impact  has  not  occurred  prior  to  bubble 
rebound.  It  further  suggests  that  the  regime  in  which  maximum  pressure  increases  as 
depth  decreases  is  characterised  by  bubble  rebound  occurring  after  the  bubble  has  evolved 
into  toroidal  form.  It  is  therefore  instructive  to  inspect  the  bubble  shapes  for  a  depth 
in  the  regime  where  the  bubble  is  toroidal  at  minimum  volume  but  the  peak  pressure  is 
decreasing  with  decreasing  depth  in  order  to  identify  any  possible  behaviour  that  may 
characterise  this  regime.  The  motion  of  the  toroidal  bubble  at  a  depth  of  20  m  is  shown  in 
figure  11.  After  jet  impact  the  bubble  continues  to  collapse  and  a  depression  in  the  bubble 
shape  is  formed  and  swept  around  the  bubble  by  the  circulation.  As  the  bubble  begins 
to  re-expand  the  depression  vanishes.  It  might  be  supposed  that  in  some  circumstances 
this  depression  would  pinch  the  bubble  off  in  the  middle  and  form  two  counter-rotating 
toroidal  bubbles.  This  pronounced  depression  may  be  a  signature  of  this  regime  where 
peak  pressure  decreases  as  the  depth  decreases.  It  is  noted  that  the  formation  of  such  deep 
depressions  in  the  bubble  shape  following  transition  to  a  toroidal  geometry  is  evident  in 
computations  of  transient  bubble  collapse  near  a  rigid  boundary  (Tong  et  aL  1999).  In 
these  cases  the  boundary  induces  asymmetry  in  the  flow  field  and  causes  jet  formation, 
rather  than  gravity  as  in  this  case. 


7  Discussion 

In  this  report  the  technique  of  Arons  (1948)  for  determining  the  parameters  that 
characterise  the  adiabatic  model  of  underwater  explosion  bubble  contents  has  been  revised 
to  take  account  of  the  effect  of  non-spherical  bubble  collapse.  The  technique  has  also  been 
amended  to  provide  increased  accuracy  with  respect  to  the  deviation  between  theoretical 
and  experimental  values  of  the  period. 

The  parameters  predicted  are  close  to  those  originally  given  by  Arons  (1948),  but 
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investigation  of  the  effect  of  non-spherical  collapse  on  the  pressure  field  suggests  that  the 
agreement  is  somewhat  fortuitous.  The  boundary  integral  method  has  shown  significant 
deviation  in  the  maximum  pressure  at  rebound  brought  about  by  translation,  non-spherical 
collapse  and  toroidal  bubble  formation.  It  has  been  shown  that  the  variation  of  the  peak 
pressure  with  depth  or  charge  mass  exhibits  a  distinct  local  minimum  at  the  point  where 
the  transition  to  a  toroidal  bubble  coincides  with  the  bubble  minimum.  It  has  further 
been  shown  that  if  the  charge  mass  and  depth  are  such  that  only  weak  oscillation  of  the 
bubble  occurs  then  jet  formation  is  suppressed  and  the  spherical  model  predictions  of  the 
pressure  are  close  to  those  provided  by  the  boundary  integral  method. 

As  it  happens,  the  charge  masses  and  depths  used  by  Arons  et  a/.  (1948)  to  compile 
their  peak  pressure  data  are  in  the  regime  where  the  bubble  is  approximately  spherical  at 
minimum  volume.  The  formation  of  liquid  jets  upon  bubble  collapse  and  the  transition  to 
a  toroidal  form  were  not  fully  understood  at  that  time.  Arons  et  al  (1948)  remark  that 
the  “twelve-pound  charge  is  the  largest  that  could  be  tolerated  at  a  distance  of  500  ft  below 
the  vessel  performing  the  experiments  An  unforeseen  consequence  of  this  practical 
experimental  limitation  is  that  Arons’  (1948)  original  analysis  is  essentially  correct,  despite 
the  more  complex  collapse  behaviour  of  the  bubble  and  pressure  field  variation  documented 
in  this  report. 

In  performing  the  work  reported  herein,  a  number  of  additional  points  have  come  to 
fight  concerning  various  bubble  parameters.  As  noted  in  the  text,  the  effect  of  deformation 
from  spherical  shape  was  determined  to  have  virtually  no  effect  on  the  bubble  period.  This 
was  established  by  performing  a  series  of  calculations  for  the  5.448  kg  charge  over  the  range 
of  depths  10-150  m.  This  is  equivalent  to  performing  calculations  over  a  substantial  region 
of  the  a,  6  parameter  space.  These  same  calculations  additionally  established  the  extent 
to  which  a  buoyant  bubble  remains  spherical  through  the  first  expansion.  The  maximum 
bubble  radius  was  computed  in  four  ways:  as  half  the  maximum  horizontal  and  vertical 
extents  of  the  bubble  and  based  on  the  maximum  volume  and  cross  sectional  area  of  the 
bubble.  The  values  determined  deviated  from  a  non-dimensional  value  of  1  by  at  most 
0.15%. 

During  early  phases  of  this  work,  attempts  were  made  to  determine  the  bubble  parame¬ 
ters  using  period  data  and  the  maximum  radius  data  of  (3).  A  typical  approach  might  be 
to  determine  7  and  e  as  done  here  and  then  fix  k  so  that  the  deviation  between  theoretical 
and  experimental  expressions  for  the  maximum  radius  is  minimised.  This  methodology 
is  not  very  satisfactory  for  a  number  of  reasons.  First,  the  theoretical  expression  for  the 
maximum  radius  shows  a  systematic  trend  with  depth  that  is  not  evident  in  the  expression 
of  (3),  after  the  dependence  upon  (£+£o)”1/3  is  accounted  for.  Second,  as  noted  in  section 
2,  there  is  disagreement  as  to  the  coefficient  J  that  appears  in  (3).  The  variation  of  this 
value  causes  large  discrepancies  in  the  values  of  k  that  are  determined  and  this  impacts 
upon  predictions  of  loading  significantly. 

Consideration  of  the  difficulties  involved  in  compiling  data  for  the  maximum  radius 
indicates  that  these  data  are  probably  not  as  reliable  for  use  in  determining  the  bubble  pa¬ 
rameters  as  period  and  peak  pressure  data.  Explosion  bubbles  are  typically  characterised 
by  a  surface  distorted  by  phase  changes  and  detonation  products.  Hence  it  is  difficult 
to  precisely  determine  the  maximum  radius.  Additionally,  it  is  difficult  to  deploy  the 
requisite  high  speed  photography  in  large  field  trials. 
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As  with  all  theoretical  studies  it  is  appropriate  that  experimental  work  be  performed  in 
order  to  assess  the  validity  of  the  conclusions  drawn  here.  There  remains  a  requirement  for 
decisive  experimental  data  concerning  the  loads  experienced  at  boundaries  due  to  nearby 
bubble  collapse  and  jet  impact.  Once  these  data  are  available  it  would  be  appropriate 
to  compare  the  results  with  computations  using  the  boundary  integral  method  and  the 
parameters  given  in  this  report.  Such  an  exercise  would  provide  valuable  evidence  with 
which  to  assess  the  practical  utility  of  the  work  reported  here. 

Final  comment  concerns  the  error  estimates  for  the  explosion  bubble  parameters,  pro¬ 
vided  here  for  the  first  time.  A  conservative  approach  to  determining  these  has  been 
adopted  for  three  reasons.  The  data  of  Arons  et  al  (1948)  show  a  reasonable  degree  of 
scatter  for  the  peak  pressure.  Second,  the  numerical  algorithm  used  for  the  non-spherical 
bubble  calculations  exhibits  a  small  degree  of  variation  as  the  various  numerical  parameters 
that  characterise  it  are  varied.  Finally,  this  work  has  neglected  fluid  compressibility.  Arons 
et  al.  (1948)  measured  the  pressure  at  around  two  maximum  radii.  It  is  not  unreasonable 
to  suppose  that  the  amplitude  of  the  peak  pressure  may  be  modified  by  compressibility 
compared  with  those  values  predicted  here  using  an  incompressible  model.  The  effect  of 
compressibility  is  to  slightly  mitigate  the  pressure  so  that  if  the  fluid  were  incompressible 
the  pressures  would  be  higher  than  those  measured.  It  is  supposed  that  the  generous  error 
bounds  used  here  for  the  peak  pressure  cover  a  sufficient  range  to  account  for  this  effect. 
The  error  bounds  provided  here  for  7,  e  and  k  are  thus  considered  reliable,  especially  in  the 
context  of  predicting  loads  to  structures  caused  by  nearby  underwater  explosion  bubble 
collapse. 
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